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SUMMARY 


The  downward  motion  of  small  paiticles  in  an  isothermal  turbulence 
free  atmosphere  has  been  studied  for  the  case  of  variabJe  dit fusion 
coefficient,  utilizing  tvjc  models.  In  the  first  of  these,  the  diffusion 
coefficient  has  been  taken  as  an  exponential  function  of  altitude  while 
in  the  Second,  the  variation  has  been  assumed  to  be  parabolic  In  the 
exponential  model  an  analytic  solution  for  simultaneous  diffusion  in  the 
vertical  and  horizontal  directions  has  been  obtained;  whi le  for  the 
parabolic  model,  an  analytic  solution  for  vortical  motion  only  has  been 
obtained.  The  ca Iculat ions, which  have  been  restricted  to  the  exponential 
case  show  the  development  of  a  profile  of  constant  shape,  all  points  of 
which  move  with  the  same  velocity,  independent  of  the  mass  cf  the  particles 
These  results,  which  arc  derivable  from  our  solution,  are  in  agreement  with 
the  numerical  work  of  Banister  and  Davis  It  is  also  shown  that  the 
qualitative  nature  of  the  results  is  independent  of  the  initial  distribution 
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CALCUIAT'ON  OF  THE  DIFFUSION  OF  SMALL  .PART  ClfS 
IN  A  NON  UNIFORM  ATMOSPHERE 


by  Milton  M  Klein  and  Kwanjj  V: 


1 .  Tnt  rodact ion 


of  cTitaiicd  1  releases  and  n;ciear 
explosions  have  stisiulaied  interest  in  the  pe’^eral  pr'^Men'  of  diff;si\e 
motion  of  particles  in  the  atmosphere-  T.ene  ra  1  analytic  scl:ticps  a  >-c 
not  available  because  of  the  variatic"  of  density  tand  he^ce  diffision 
coefficient)  in  the  atmosphere  Further  ccmplcxity  is  added  tc  the 
problem  because  of  the  general  downward  drift  d  ■e  to  the  p  ra';  i  * -a  t  icna  1 
field  If  the  region  of  interest  is  net  toe  large  the  diff.sin 
coefficient  may  be  considered  constant;  and.  fc'  iht  rase  of  negligible 
drift  velocity,  standard  metbeds  of  soiation  may  be  app  1  i cd  '  ' ^ ^ ^  The 
effect  of  a  drift  velocity  has  been  taken  into  acc;:nt  for  the  case  cf 
constant  diffusion  and  ccnslani  drift  velocity  by  Cliandr-isc  kh  a  r  wbc  se.? 
a  change  of  variable  tecbniqje  tc  reduce  I  be  ..-r.  blcm  t  a  sla^'dard  c-'e 
in  heat  conduction 


An  invest  igat  icn  of  the  desc<  nt  rf  small  p-irt,r]es  th-n  gb  a- 
exponential  atmosphere  t-y  r'cars  cf  n' .m':  r  i '  i  i  sol  'i-n  pf  t  bo  gr'e'sint' 
equations  has  recently  f^rn  made  by  banislc-r  -and  r;j.,  1  s  the  st  dv 

■••as  icbiticifcd  i  j  mei  I  cn  in  the  vertical  plar-e  j^e  the  st-iking 
results  of  tbeir  work  is  Iht  change  in  a  sbert  time  inte'val  f  t hr 
initial  density  distribution  into  a  censtant  pr  file  as  i'  desttods 
through  the  atmosphere  An  analytic  scliticn,  metivated  bv  this  resilt 
has  been  developed  by  t-ranri^w i'l  the  form  of  an  infinite  sc’^ies  ct 


associated  Lagueri'  polynomials  whose  coefficients  are  obtainable  from 
integrals  over  the  initial  distribution 

Because  of  the  general  interest  in  the  problem  of  atmospheric 
diffusion  and  the  scarcity  of  detailed  solutions  when  the  non  uniformity 
of  the  atmosphere  is  taken  into  account,  we  have  studied  the  problem  in 
detail  for  two  models.  In  one  of  these  we  have  taken  the  diffusion 
coefficient  as  an  exponential  function  of  altitude  'while  rn  the  other 
one  the  variation  is  parabolic.  For  the  sake  of  simplicity,  the  atmos¬ 
phere  has  been  assumed  isothermal.  Ihe  exponential  model  is  chosen  in 
accordance  with  the  approximate  exponential  variation  of  density  of 
the  atmosphere  with  altitude  and  is  convenient  for  U'S themat icai  analysis 
when  considering  an  infinite  or  semi  infinite  region  The  parabolic 
model,  while  it  can  be  used  for  a  semi-infinite  space,  is  more  useful 
for  a  finite  region  where  the  analysis  for  the  exponential  model  becomes 
awkward  and  inconvenient 

We  shall  restrict  ourselves  in  the  present  analysis  to  the  case  of 
molecular  diffusion  and  therefore  will  not  consider  in  any  detail  the 
effect  of  turbulence-  Since  the  atmosphere  is  reasonably  turbulence 
free  down  to  altitudes  of  about  100  km,  we  may  assume  our  calculations 
are  applicable  down  to  100  km  Any  further  extension  be  .ow  this  altitude 
should  be  examined  with  extreme  care,  The  effect  of  a  turbulent  layer 
upon  the  particle  distribution  may  be  analyzed  Dy  methods  similar  to  the 
one  presented  in  Refers  .ce  4, 
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2 .  Analysis  and  Formulation  of  Model 


We  consider  a  non-uniform  atmosphere  in  which  the  diffusion 
coefficient  D  is  a  known  function  of  altitude  z  A  foreign  gas 
deposited  in  the  atmosphere  will  spread  in  all  directions  due  to 
diffusion  and  vjill  drift  downward  under  the  action  of  the  gravitational 
field-  The  drift  velocity  v.'ill  be  controlled  hy  the  frequency  with 
which  the  gas  molecules  collide  with  the  atmosphere  The  basic 
diffusion  equation  and  the  dependence  of  the  drift  velocity  upon  the 
collision  frequency  may  be  obtained  from  the  Boltzmann  Equation 
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the  internal  collisions  among  the  gas  particles  while  I 


6t/  ext 


the  effect  due  to  collisions  between  the  gas  and  ambient  air  particles 
Since  the  gas  particle  density  becomes  quite  low  after  a  short  period 

If 


of  time,  we  shall  neglect  the  self-collision  term  I  tt  I  ■  .  Ir. 
our  case,  the  only  external  force  is  the  gravitational  field  which  acts 
in  the  negative  z  direction  so  that  Equation  (1)  may  be  written  in  the 
form 

li  . 


c)f  ,  ■  Df 

dt  - 

dr 


-1 


(2) 


wh 


ere  c^  is  the  .  'omponent  of  the  velocity 


4 


We  now  obtain  our  macroscopic  equations  of  motion  by  multiplying 
Equation  (2)  by  an  appropriate  power  of  the  velocity  and  integrating 
over  the  velocity  space.  For  our  purposes  it  will,  be  sufficient  to  take 
the  zeroth  and  first  powers  of  the  velocity  We  then  obtain 
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dc 
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f  '/ffN 

d  c  -  g/  ^  c  d  c  =  ^  c  c  c  (4) 

dr  y  J  3  -  '  " 


In  Equation  (3),  since  integration  over  velocity  is  independent  of  time 
and  space,  we  may  write 


a  “■  3  / 


-  b'r; 


at  '  ^  ^  ^  =  dt 


where  the  particle  density and  the  macroscopic  velocity 
by 


are  defined 


( 5i 


Since  the  velocity  distribution  vanishes  at  the  limits,  tlic  term  in  — - 
gives  zero  contribution.  The  collision  integral  in  Equation  (3)  may  be 
shown  to  vanish  when  integrated  over  the  velocity  space  The  macroscopic 
form  of  Equation  (3)  therefore  gives  the  usual  equation  of  continuity 

7  •  (flv)  ==0  (61 
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it  is  convenient  to  consider  s.  typical  compouenc 
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For  ■  =  3  we  get 


and 
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!  ^ —  c.  d  c,  - 
!  dc -  3  3 

J  J 


Sf 


.1 


■-'3  ^ 


-  Jf  d  c. 


d  c  =  -  Jf  d  c  = 

We  may,  therefore,  write  Equation  (4)  in  the  form 


(7) 


r- 

V  ■  c  f  d  c  +  g 
C  C 


*.3  =/(!{) 


ext.  i  (8) 


Because  of  the  complicated  form  of  Equation  (6),  the  direct  solution  of 
our  problem  by  means  of  Equations  (6)  and  (8)  is  a  formidable  task.  We 
shall,  therefore,  utilize  Equation  (8)  to  give  us  an  approximate  expression 
for  the  drift  velocity  and  then  solve  Equation  (6).  In  Equation  (8)  we 
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note  that 


-T—  I  is  prooortioned  to  f  and  to  the  collision  frequency 

ct  J  ext  - 

which  is  a  function  of  a  velocity  and  position  The  dependence  upon 
the  velocity  distribution  is  quite  weak^  and  may  be  ignored  here;  we 
may,  therefore,  write 


r/fiA  ^  r  r  .  - 

;  fr— 1  c  d  c  ^  /  V  (r,  V)  f  c  u  c 

X^tjext.  a  J  ^ 


=  4’, 


f  c  d  c 


where  v  is  the  average  value  of  the  collision  frequency  over  the  velocity 

space  The  velocity  v^  is  the  sum  of  the  diffusion  velocity  v^^  and  the 

drift  velocity  v  v  j  ■  i  c  n 

cyg.  Since,  the  drift  velocity  goes  to  zero  for  g  =  0, 

We  may  identify  the  drift  velocity  term  on  the  right-hand  side  of 

Equation  (8)  with  the  term  containing  g  on  the  left  side  We  therefore 

have 


og 


(9) 


from  which 


\X  J 

V  -  - 

Q-g  V 


(10) 


We  now  express  the  velocity  v  in  Equation  f6)  es  the  sum  of  the 

diffusion  velocity  v„  and  the  drift  velocity  v  and  write 
D  •'  g 


(11) 


/ 


where,  for  our  case,  v  has  only  a  z  component  v-  ,  given  by  Equation  (10) 

g  Jg 

The  diffusion  velocity  Vj^  is  given  by  Fick's  law^^^ 


Vi  = 


(12) 


where  D  is  the  diffusion  coefficient.  We  may,  therefore,  ivrite 
Equation  (11)  in  the  form 


^  -  i'  •  (D  V^/)  +  V  ■ 


)  =  0  (13) 
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Analysis  of  the  Diffusion  Equation 


The  solution  of  Equation  (13)  for  the  general  three-dimensional 

case  is  quite  difficult.  It  is  possible  to  s^a rate  Equation  (13)  into 

two  equations;  one  governing  the  motion  in  the  horizontal  plane  and 

the  other  governing  the  motion  in  the  vertical  plane,  provided  we  assume 

that  the  coupling  between  the  tiv’o  motions  is  weak  The  significance  of 

this  approximation  and  an  interpretation  of  the  resulting  solutions  will 

be  presented  later.  A  convenient  and  useful  variable  whose  differential 

equation  describing  the  vertical  motion  is  not  too  complicated  is 

provided  by  the  particle  density  integrated  with  respect  to  x  and  y, 

i.e.,  the  number  of  particles  per  unit  z  length  To  obtain  this  equation 

we  first  express  Equation  (13)  in  cartesian  form  and,  noting  that  v  is 

8 

a  function  of  z  only,  obtain 

^  ■  fc  <“  -  ly  <■>  l?>  -  fe  <»  0  ■ 

Integrating  with  respect  to  x  and  y  from  to  -l-«  and  noting  that 
and  vanish  at  the  limits,  we  obtain 


5n  dn 

■VT~  -  ^  (D  —  -  v  n  )  =  0 
ct  dz  dz  g  V 


.  ^  r  f-h.  ^ 

"1  “  J.c= 

is  the  numbei  of  particles  per  unit  z  length.  We  now  look  for  a  solutic 
where  n^  (z)  describes  the  vertical  motion  and  write 
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‘ij(x,y,z')  -  nj(z)  n2(x,y,z) 


(  16) 


where,  since  we  are  assuming  weak  coupling,  has  only  a  weak  z 
dependence-  Substituting  Equation  (16)  in  Equation  (13)  yields 


n, 9n_  n~3n 

i-  — —  =  n  "  (D  V  n  -  V  n  )  +  V  (D  n,  7  n,) 

3t  3t  2  1  g  1  12 


+Vn-  (OPn  vn) 

2  1  g  1 

which,  by  virtue  of  Equation  (15),  can  be  simplified  to 


117) 


•‘1"“2 
at 


-  =  V  ■  (D  n.,  V  n„)  +  17  n„  (D  7  n  v  n  ).  (18) 

12  2  1  g  1 


Since  D  and  are  functions  of  z  only,  wo  may  write  Equation  (18)  in  the 
scaler  form 


n,an„  ^  an,  an  dn,  /  3  n„ 

-ir^  ■  17  -57>  *  '■’ir  -  “g".'  57^  *  ^ 


In  view  of  the  assumption  of  weak  coupling,  we  may  neglect  the  terms 
involving  the  derivatives  of  n2  with  respect  to  z  to  obtain 


at 


2  =  0 


/  ^""2 


(20) 


(19) 


\  ax  cy  / 

Since  D  is  a  function  of  z  only,  Equation  (20)  represents  the  diffusion 
of  material  in  the  x-y  plane  for  a  constant  diffusion  coefficient  whose 
value  is  determined  by  the  particular  altitude  at  which  the  diffusion 
is  taking  place.  A  solution  of  Equation  (20)  for  the  case  of  a  source 
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in  the  xy  plane  is  given  by 


(8) 


"2  = 


r  2  2 

4Dt 


(21) 


where 


CC  ''C 

„=/  / .. 


n^  dxri  S' 


is  the  source  strength  Because  of  the  drfinition  of  n^  (Equation  (lo) i. 
the  value  of  K  for  our  case  equals  unity 

At  thrs  point,  we  consider  the  significance  of  the  assuirpt  ion  of 
•.;cak  coupling  If  the  motion  is  uncoupled,  tht  n  in  a  time  t  a  particle  will 
diffuse  a  given  amount  in  the  x  y  plane  independent  of  its  aownvv’ard 
motion  Actually  the  particles  move  in  a  curbed  path  toward  lower  values 
of  the  diffusion  coefficient  so  that  the  extent  of  travel  in  the  ;c-y 
plane  is  actually  less  than  that  obtained  under  the  assumption  of  no 
coupling.  Thus,  the  diffusion  front  calculated  under  out  assumption  of 
weak  coupling  will  have  a  larger  extent  in  the  x  y  diiection  than  the 
true  diffusion  front  The  extent  of  diffusion  along  the  z  axis  will  be 
correct  since  Equation  (15)  for  this  motion  is  exact 

The  three  dimensiorial  diffusion  equation  may  be  treated  from  an 
alternative  point  of  view  We  first  wi ite  Equation  (13)  in  the  form 


ill 

3t^ 


A  ;7V; 

~  (D 
Sz  uz 


I 

1" 


V  7/  )  D  -  +  --- 
8  ' 


--j-  =0  (22) 
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ov 
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and  then  Fourier  transform  Equation  (22)  with  respect  to  x  and  v  to  obtain 
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Sn  ^  an 


V  n^)  ■  D  =  0 

g  f  -■  ■  f 


where 


r 

r "" 

/  ***  r 

/  ^ 

(r  ,t)  exp|  i(x5  : 

yd)j 

•"  ■’  rr,  0 

and  wc  have  made  use  of  the  boundary  condition  that  ,  dx  ,  and  Sy 

vanish  at  infinity  We  note  that  for  the  special  case  5  =  =  0>  ihe 

transformed  density  function  n^  reduces  to  the  density  along  the  z  axis  , 
n^Cz)  Ksnce,  the  density  n^fz)  is  obtainable  as  a  special  case  of  the 
Fourier  transform  function  n^  We  shall  work  with  Equation  (23)  rather 
than  Equation  (15)  and  see  to  what  extent  we  can  solve  for  the  general 
three-dimensional  density"^  without  appealing  to  separation  of  variables 
and  a  perturbation  procedure  for  obtaining  the^  diffusion  in  the  horizontal 
plane, 


At  this  point,  it  is  convenient  to  express  our  equations  in  dimension 
less  form;  we  therefore  set 


The  Fourier  transform  of  Equation  (26)  with  respect  to  x'  and  y 


yields  in  place  of  Equation  (23) 


T  d_  4  ^■^f 

t- '  '  :>  i-,  '  ( 


at'  ^2  Sz 


az  ■ 


J-  Vf)  ^  ^  "f  "  ° 


C2i 


where  n^  is  given  bv 
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cc  00 

=  f  f  ‘’J](r,t’)  exp  [i(x's'  +  y  ' "  ‘  )]dx V ' 

■^20  *  OO 


(28) 


We  shall  take  the  z  axis  as  positive  upward  so  that  the  drift  ve I 


in  Equation  (27)  is  negative  We  therefore  replace  v^  by  u^  where 

u  is  positive.  ft  is  also  convenrent  to  express  u  in  terms  of  D  by 
g  g 

noting  chat,  from  elementary  kinetic  theory 

2 


^  1  ,  1  V 

D  =  3  V  ,.  =  3  - 


(29) 


where  v  is  the  mean  thermal  speed  and  the  mean  free  path  Since  u 
from  Equation  (10)  =  g/p,  it  follows  that 


g 


D  2 


(30) 


and  is  constant  for  an  isothermal  atinosphere 


We  may  then  write  Equation  t27)  in  the  form 


Bn, 

_ I  r 

ot ' 
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2  Bz' 


/  An 


DI  ,.,2 
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,.2 


)  =  0  (31) 
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4,  Exponential  Distribution 


We  consider  now  the  analysis  of  Equation  (31)  for  the  case  of  an 
exponential  distribution  of  diffusion  coefficient  with  altitude  We 
take  the  origin  at  the  earth's  surface  and  unite 


f.e-'”  ,32, 

o 

where  D  is  a  reference  value  at  the  origin  and  II  is  the  scale  height 
o 

Since  we  have  assumed  an  isothermtil  atmosphere,  H  will  be  taken  as 
2-23nct-pr.t  Tho  differential  Equation  (31)  then  becomes 


f  e  \  1  '' 

~  -r  D  —  exp  (z/H)l  ^  L  ^  n  J  -s 
Ot  ,  /  o  \cz  3  f /  1  r 

L  \  V  y  ,1  L 


2  2 

exp  (z/H)(y  T,'  )  =  0 


We  note  that  Equation  (33)  simplifies  considerably  with  regard  to  the 

coefficients  if  vjc  take  If  as  the  unit  of  length  1  and  H^/D  as  the  unit 

o 

of  time  T;  v.'g  then  have 


oil-  I  CU-  - 

ne-T  •  exp  (z‘)  I  --T  +  H  -§  n 
ct  cz  \dz  2  t 

1  '  V 


exp  (z  ’)  ( bj  =  0 


2  3kT  kT 

Noting  that  v  - - and  11  = - where  m  and  m  are  the  masses  of  the 

m  m  g  a 

d 

foreign  gas  and  the  atmospheric  molecules;  and  dropping  primes  for 
convenience  for  the  remainder  of  Section  4,  we  may  write  Equation  (34)  as 


3t  dz 


9  f  f'" 

--  exp  (z)l  ■; 


—  +  an 
z  f, 


2  2 

exp  (z)(c  -f  T]  )  n^  =  0  (35) 


V-  gH  ^  , 

where  a  =  measures  the  ratio  or  the  masses 

2  v  D 


I 


We  now  take  the  Laplace  transform  of  Equation  (35)  and  obtain 


which  reduces  to 


^  *  (a  -  n  —  ’  r  a  -  -  p  exp  ( -z)!  M  =-t(z 

dz^  clz  L  J 


z  )  fAO' 
o 


by  use  of 


[-> 


M  -  --  N  exp  -  i  fx  V  -  y  ) 


oj 


We  now  change  the  independent  variable  z  by  J  -  exp  -z/2 ■  to  obtain 

2  ,2 


r2  cl^M  1  .  dM  .  ,  , 

C  — :7  -  <2a  '  I )  77-  •  4(a 

dC^ 


-  pC")M  ^  -2^  6t£  -  '  ' 

o  o 


'41 


and  further  introduce 


.-(a-  1  ). 


to  yield 


.2  O  r  ^ 

"  ..2  '  dC 


<5^  -  'ip  "  •2£  6'C  -  £  I 

'o  o 
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where 


2  2  2  2 

3  -  (a  •  1)  >  4(:  •  P  > 


Since  Q  -  exp  (  -z/2)-  Q  is  defined  for  0  <  £  <  1  We  now  define 
the  analytic  continuation  of  0  into  the  region  £">  1  by  Equation  i'42; 

The  right-iiand  side  of  Equation  (42'>  is  zero  at  £  -•  I  so  that  we  shall 
require  Q,  its  first  and  second  derivatives  are  continuous  at  £  -  1 
and  Q  —  0  as  £  —  ~  This  corresponds  to  the  assumption  tiiat  the 
Dc.rticle  density  approacties  zero  as  z  —  •  Although  the  pliysical 

spare  covers  only  positive  values  of  z,  the  pai'^irle  density  is  not 
significantly  affected  by  the  boundary  conditions  at  z  -  0  until  the 
majority  of  particles  itacii  iiie  ground  level  Furthermore,  the  present  e 
of  turbulence  below  100  km  renders  invalid  any  molecular  diffusion  model 


belovj  this  level  and,  in  addition,  the  turbulent  layer  behaves  like  a 
particle  sink.  We,  therefore,  expect  our  present  solution  to  be  a  good 
approximation  to  the  actual  situation  down  to  the  iOO  km  level. 


in  the  p -plane  is  then 


N(C>P)  =  “3  exp  I  i(x  §  i-  y  r 

H 


a  T-i  •j^{2i\fpQ)  H,3^^^(2i'/pr 


)  ecr 

-o  ~o 


<-  Jq(2i''pr  )  Hq^'\2iVpr)  aCC  ■  ■'  )]  1^6) 


Applying  the  inverse  Laplace  transform  to  M(c,P) 

,  .  ‘  n 


2-i  >io.,c 


N(»,P)  exp  (pt)  dp,  c  ^  o 


(47) 


we  obtain 


.  ^2  2 


=  —  exp  ygh) 


1  •  a  i-i-a 


exp  I  ' 


(48) 


The  particle  density^  (r,t)  is  then  given  by 

(r,t)  ==  ®-'^P  [  ■*■  yT>j  ded- 


CO  00 


4.2  , 


I  dc  exp  i  iCx  -  '^^P  [•i(y  *  yQ^n] 


M  1  ..  :  r  2  ]  ^2r;  \ 

o  I'a  .i-:-a  ,-o  I,  ='’o 

.X'Co  c  exp^ — r — 


Setting  X  -  X  =  r  cos  6,  y  -  y  -  r  sin  0 
o  o 


T*  =  p  sin 


(50) 


and  using  the  integral  representation  of  we  have 

j  dg  exp  ^-i(x  -  j  ‘I.  exp|-i(y  -  y^)pj  I 


2  ?  7 

45  ■  ‘‘T  ■  ia  •  I  I 


a 


(  '  >  J  (  cT  •  ci  r 
o 
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2rr 

where  j  =  and  a  =  a  -  1.  The  particle  density  may  therefore 

be  written  as 
n 


,2 


! 


7|(r,t)  2  3 

4^  H  t 


1-a  ,1+a  o  I  ; 

»  e^-Pl  - - 1 - ■  ^  ,  2  2 

^c,  .4c  -Q 


o 


/^aV( 

\  I' 


(x  ■  'a)‘ 


(y 


.  2 


(5  1) 


The  integral  occurring  in  Equation  (51)  is  a  Hankel  transform  of 
order  zero  and  can  ea.sfly  be  evaluated  in  anv  particular  case  As 
previously  indicated,  the  onc-d imcns iona 1  solution  n|(z,t)  may  be  obtained 
by  setting  .  =  .  =  o  in  the  solution  for  n,  We  then  obtain  from 

Equation  (48) 


*■>  .  I  . 

n  j  a  -  l-:-a 

=  iiT 


',2 


'  I  1 

/  a  ■  1 


'iiioN 


(52) 


Some  useful  information  concerning  the  late  time  history  of  the 
particles  may  be  obtained  from  the  solution  for  njfz.c)  Tor  large 
values  of  t,  Che  Bessel  function  in  Equation  (52;  may  be  expanded  in  a 
power  series;  and,  upon  retention  of  the  leading  order  term,  we  obtain  as 
the  asjmipcotic  solution  for  n|(z,t) 


n  /.2 

,  o  1 

"l^^'  ^  "  H  r(a)  V  t 


(0  ... 


(52a) 


The  exponential  function  is  ciose  to  unity  for  large  t  so  that  dependence 

of  the  density  function  n  upon  z  and  t  occurs  only  through  the  sineie 
2  ‘  ' 

variable  £  .  The  altitude  z  is  related  to  the  lime  t  by 
t 


=  exp  (  -z)  =  ct 


(52b) 
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I 


where  c  is  a  constant  which  depends  upon  the  density  nj  and 
a,  H,  and  2^.  We  shall  find  Equations  (52a)  and  (52b)  usef 
the  results  of  the  calculations. 


parameter  s 

il  in  interpret  ing 
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5 .  Parabolic  Model 


We  examine  now  the  case  where  the  diffusion  coefficient  is  a  parabolic 
function  of  altitude;  it  is  convenient  to  write  the  diffusion  coefficient 
in  the  form 


rv  ^  2  \“ 

D. 

1 


(53) 


where  b  is  the  altitude  at  which  the  gas  is  released  and  D,  is  the  reference 

J. 

value  of  o  at  that  altitude.  Substituting  Equation  (53)  in  Equation  (31) 
yields 


I  !i 


L"  b" 


T 


(5i) 


Parallel  to  the  analysis  for  the  exponential  case,  we  now  choose  b  for 
the  unit  of  length  L  and  “  as  the  unit  of  time  T;  we  may  then  wite 
Equation  (54)  in  the  form 


5t' 


( - 


,2 


n^  -  0  (55) 


If  we  apply  the  Laplace  transformation  to. Equation  (55),  the 

inhomogeneous  term  of  the  resulting  ordinary  differential  equation  is  a 

delta  function.  Because  of  some  difficulty  in  finding  the  corresponding 

Green's  fun  cion,  we  have  explored  the  use  of  the  Lagrange  variable 

z  ' 

hf  =  /  n^dz'  (56) 

-'o 

for  which  the  inhomogeneous  term  is  a  step  function  Because  of  the 
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presence  of  the  c'  and  t'  terms  in  Equation  (55),  however,  the  transformed 
equation  i c  an  integral  rather  than  differential  equation  We  have  there¬ 
fore  restricted  ourselves  in  the  present  analysis  to  the  density  function 
nj^(z,t)  along  the  z  axis.  Setting  =  "i'  =  0  in  Equation  (54)  and 
integrate  vjith  respect  to  z  ‘ ,  we  obtain 


dh 

dt 


r  -  (  z '  —~2  ^ 


^  bz'" 


0 


(57) 


V 


where 


c  . 

h  =  J  nj^(z  , 


'  .t)dz  ' 


(58) 


and  we  have  imposed  the  boundary  condit ion  that  the  flux  vanishes  at 
z  =  0 


Applying  the  Laplace  transform  to  Equation  (57)  and  omitting  primes 
for  the  remainder  of  Section  5  yields 


2  ^  2  ^ 


dz 


2  ■  2 


pH  -  -h(z ,0) 


(59) 


where 


"  =  /' 


h(z,t)  exp  (-pt)dt 


(60) 


and  h(z,o)  is  the  initial  distribution  in  terms  of  h  An  exact  solution 
of  Equation  (59)  Leads  to  hypergeometr ir  functions  which  introduce  some 
difficulties  in  obtaining  the  Laplace  inversion  We  have,  therefore, 
simplified  Equation  (59)  by  approximating  z  in  the  —  (drift  velocity) 
term  by  cvz  where  a  is  a  suitable  avei age  of  z  over  the  interval.  Since 


we  are  interested  in  values  of  the  altitude  above  100  km,  we  expect  the 
approximation  to  be  a  useful  one.  Wo  may  then  write  Equation  (59)  as 


where 


3  L. 


2  a  H  ^  B  MI  „  ,  ,  , 

z  — ~  -  pH  =  h(z,o) 

dz“ 


Solving  Equation  (61)  for  the  complementary  function  H 


we  obtain 


m  m., 

H  =  Az  Bz 

o 


whore  m^  and  m^  are  given  by 


mj  =  c  ,  V  c  +  p 

m2-c  '/c^  +  p  (63) 

and  c  =  ^""2'-"  .  For  the  initial  distribution,  we  choose  a  delta  function 
n 

at  z  =  1,  — ^  f(z  -  1)  The  corresponding  Lagrange  function  h  is  then  the 
Heaviside  step  function  8(1).  A  particular  solution  H^  of  Equation  (61) 
for  this  inhomogeneous  term  is 


Hp=^?(l)  (6A) 

and  the  complete  solution  H  =  H  +  H  is  then 

o  p 

m  m2 

H  =  Az  +  Bz  ,  z  <  1 
m  m„  n 

H  =  Cz  -1  Dz  z  >  1  (65) 

where  we  have  used  different  constants  in  the  tvro  regions  because  of  the 
step  function  at  z  =  1. 


We  consider  now  the  boundary  conditions  for  our  problem.  We  shall 


I 


place  perfect  reflectors  at  z  -  0  and  at  an  upper  level  z  =  d  and  there¬ 
fore  require  zero  flux  at  these  locations.  Froni  Equations  ''61)  and  (64) 
the  flux  F  is  given  by 


t  =  const,  (n.i  n  (  1) 
•  p  o 


In  addition,  wo  shall  make  use  of  the  continuity  of  n.|(z,t)  and  hence  of 
h(z,t)  anc:  H(z,p)  for  t  ^  0  at  z  =  1,  Since  the  real  part  of  the  parameter 
P  is  positive,  the  exponent  m2  in  Equation  (65)  is  negative  Ho  therefore 
require  the  coefficient  B  in  Equation  (65)  to  vanish  to  avoid  an  infinite 
solution  at  z  =  G.  Application  of  our  boundary  and  continuity  conditions 


then  yields 


H  =  —  c  1-  q  -  d^~’(c  -  q)  ^c-:-q  ^  ^ 

"  I’  2qd2^ 


P  P 


z  '  1 


,  2  2 

where  q  =  c  +  p 


The  inverse  Laplace  transform  of  H  may  then  be  obtained  from  a 

( 9) 

standard  table  of  inverse  Laplace  transforms.  Several  of  the  inverse 

transforms  are  expressed  as  Faltung  Integrals  which  can  be  related  to 


Similar  relations  involving  higher  powers  or  may  be  obtained  from 

^'2. 

Equation  (68)  by  differentiating  with  respect  to  b  The  particJe  density 
nj(z,t)  is  then  obtained  from  h(z,t)  by 

.  .  <^h 

II  rz,t;  =  — 

1  oz 

The  results  obtained  may  be  expressed  as  follows: 

^C\'r  •’  1 


n  4  \  '  'Z 


.1  c.) 


(69) 


1  =  V  I2  -  -!• 


I2  =  '  ryi:tfc  (r,^  -  a,)  -  i— ;  Erfc  (a^  a^) 


-  fl)‘ 


/j2  yc 

(.,0  ai)  /  Erfc  (a 


o  ^1^ 


I3  =  —  Erfc  (a2  -  ay  -  — •  Erfc  (32  +  ay 


-  “  Erfc  (a2  +  a^^)  +  — ^  Erfc  (a 


2  - 


5 

8z 

02 


2\c 


/d^y  1 


f  .  '2 

-(a^  I-  a^; 


O  1 


-c  El  fc  (a  -f  a  )  T-  c  Erfc  (a  -  a  ) 
o  1  o  1 


^  ^  (^2  "T-  Erfc  (32  -  a  ) 

2  7. 
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A  11 
^  Sz  “  '  c 

2  V  '  t 


Ca^  -  aj) 


-  ( Z  7  4-  ^  ^  \ 

'  2  1^  :  ^  —  exp  - 

'  ■  z  ‘^'nt  I 

Er/c  (a^  +  a^)  -1-  -^-  Erfc  (a^  - 


a  2  I  d" 

o  /it  z 


a^  =  ct 


2  1.1 


2  ~  4t 

for  z  >  1, 


-.-n  — 
z 


c-  1 

n  44 
o 


,  c>J  .  , 


'•■■'here 


■'  ■  -'i  *  -■2  ■  J3  -  J4  ;  J,  -  ',,  J2  -  i,  . 


(70) 


J3  z  Erfc  (a  '  -  a  )  -  Erfc  (a,'  +  a  ) 
z  1  2  r 


J,  -  z  Erfc  (a 


2'  -;-  ap  +  z  Erfc  (a^ '  - 


ai 

dz 

^  bl 

ai2 

az  ~  " 

Bz 

] 

52 

—  Q 
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^7. 


■(^2'  ■■■ 


cl  r  ,  2, 

exp  I  -(32  •  a^l 

z  /'’t  I  J 


+  -  Erfc  (a  -:-  a  )  -  Krfc  (a  a  ) 

-  c  2  i  -  c  2  1 

z  z 


.2  1  , 

a^  =  —  -i-n  z 

2  'it- 


Although  the  solution,  for  the  par.'>ho(ir-  case  has 
analytic  form,  the  asymiptotic  solution  for  large  t  is  sufficiently  simple 
to  draw  some  interesting  and  useful  conclusions,  The  rorm  of  the  solution 
depends  critically  upon  the  value  of  c;  for  c  positive  the  leading  order 
term  depends  upon  z  but  not  t;  while  for  c  negative,  the  time  independent 
leading  order  term  cancels  out  and  the  solution  becomes  time  dependent 
The  results  obtained  for  these  two  cases  may  be  expressed  in  the  form 


_1  _  2c  2c  - 1 
1  ~  ,2c  ^ 


,  c  >  0 


(71) 


^  ^  X  X 

—  =  const  — r  exp  (ct) ,  c  =  0 
n  -i 

o  zt 

For  c  >  0  the  solution  is  in  the  form  of  a  power  law  distribution  of  the 
particle  density  with  respect  to  altitude.  This  result  is  due  to  the 
approximation  in  the  drift  velocity  term;  if  the  approximation  had  not 
been  made,  we  would  have  obtained  an  exponential  behaviour  similar  to 
that  obtained  in  Reference  3.  This  is  not  a  serious  error,  however,  since 
we  are  not  interested  in  the  region  z  close  to  zero.  We  may  see  the 
effect  of  our  approximation  more  closely  by  directly  solving  the 
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differential  equation  for  large  t.  Noting  that  the  time  derivative 
becomes  negligible  for  large  t,  we  may  write  Equation  (54)  in  the  form 


3  /  2  S"l  3gb  2  \  - 


which,  upon  noting  that  the  flux  vanishes  for  large  L,  may  be  integrated 


2  -""l  3gb  2 
^  3z  ■■■  2 


If  we  approximate  z  in  the  drift  velocity  term  by  az,  we  obtain 
2  <^"1  „ 


"i  = 


Using  n^  = 


n^dz  to  evaluate  the  constant  A,  we  obtain 


I  -  c 

A  =  •  n 

dl-i  o 


so  that 


"1  1  -  S  -S  2c  2cl 

—  =  -7-^  z  =  -n;-  z  (75) 

o  d  d 

in  agreement  with  Equation  (71) .  If  we  do  not  make  any  approximation 
in  the  drift  velocity  term,  however,  we  obtain 


2  "“1  3cb  2 

"  dr  “  ■  'r  ^  "1 

V 


=  const  =  exp  ^"2 

L  V 


which  is  an  exponential  distribution  similar  to  that  obtained  in 

Reference  3  for  the  case  of  constant  D  and  constant  u  . 

8 

The  case  of  no  gravitational  field  (3  =  0,  c  -  1/2)  yields  a 
particularly  simple  result  Setting  g  =  0  in  Equation  (73)  or  directly 
from  Equation  (75),  uc  obtain,  as  anticipated,  the  uniform  distribution 


^  =  i 

n  d 
o 


(77) 


The  case  c  =  0  is  not  directly  obtainable  from  Equation  (73)  since 
the  dropping  of  the  time-  derivative  is  not  valid  here  in  the  neighborhood 
of  the  origin  This  behaviour  may  be  seen  clearly  from  Equation  (71) 
where  we  note  that  n.^  -  0  for  large  t  except  near  z  -  Q  We  thus  have  a 
piling  up  of  particles  at  z  =  0;and  since  the  total  number  is  constant, 
a  delta  function  distribution  is  obtained  for  t  =  -'  This  result  is 
not  of  great  significance  since  it  is  due  to  our  approximation  for  the 
drift  velocity  term.  It  docs  show,  however,  that  in  a  situation  where 
the  gravitational  field  becomes  significantly  larger  than  the  diffusion 
effect,  such  a  piling  up  of  material  may  be  expected 
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6 .  Results 

Because  of  the  greater  simplicity  of  results  and  because  of  norc- 

isnmediate  applicability  to  the  physical  atmosphere,  wc  have  confined  our 

numerical  calculations  to  the  exponential  model  In  addition  we  have 

restricted  ourselves  in  the  present  investigation  to  the  one  dimensional 

density  function  nj^(z,t)  Numerical  calculations  for  the  three  dimci: 

sional  density  function  (r,L)  will  be  presented  in  a  futuie  investiga 

tion.  The  density  function  n  (7.,t')  is  prcsfnted  as  a  function  of  alii 

D 

tude  z/H  for  several  values  o£  the  time  --  t  and  the  mass  ratio  a  in 

if 

Figures  1  through  4.  Tire  scale  heiglit  H  was  taken  a-  I’O  kn  and  altitude 
z 

of  release  —  as  ?0  Tlie  figures  show  that  the  density  profile  becomes 
bi 

increasingly  narrow  as  a  increases,  and  the  rnaximuin  value  of  n^  increases 
with  a  These  results  may  be  anticipated  in  view  of  tlie  stialler  momentum 
changes  and  consequent  lesser  dispersion  of  the  heavier  particles  the 
time  taken  for  a  given  profile  to  desceno  io  a  given  altitude  is  easily 
obtained  from  the  figures.  As  an  example,  for  a  =  1  iFiguie  I),  it  take.- 

about  10  units  of  —  t  for  an  injection  at  400  km  to  descend  to  the 

«  2 
100  km  level  For  a  representative  value  of  -  400  cm  /sec,  the 

g 

corresponding  ph'/sical  time  is  10  sec  or  about  three  year'  a  reason¬ 
able  result. 

The  most  striking  feature  of  the  density  profile  is  its  shiftinc  by 
a  constant  amount  as  the  time  parameter  changes,  the  piofile  remaining 
unchanged  in  form.  It  thus  appears  thsit,  aftei  initiai  time  pf r roa 
required  to  establish  the  profile,  all  particles  mo've  witli  the  same 
velocity,  Tnese  results  have  been  obtained  previously  in  the  uumpricai 
Cax cu r a c r C'i'i s  of  Sanxster  and  Davis 
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The  development  of  a  constant  profile  and  its  consequences  are 
easily  obtcdnable  from  the  relation  between  r  and  t  given  by  Equal. ion  (52b; 
Fora  fixed  value  of  c,  we  may  write  this  relation  in  the  form 


z  I  -  Z2  =  -£n  -- 


where  2^  and  z^  are  typical  altitudes  and  t^  and  t^  the  corresponding 

times.  For  our  calculations,  the  ratio  ot  two  consecutive  times  7 — 

^1 

has  been  chosen  as  100  .so  that  the  corresponding  difference  Zj  •  z^ 

=  '-n  100  =  46  A  check  of  Figures  1  through  4  shows  that  the  altitude 
difference  between  corresponding  points  on  adjacent  profiles  i  .s  always 
about  4.6  units  The  time  required  to  establish  a  constant  profile  i.s 
easily  obtained  from  Equation  (52).  Since  the  aigument  ~~ —  ot  this 
Bessel  function  must  be  small  enough  to  yield  the  constant  profile 
solution  given  by  Equation  (52a),  the  required  time  must  be  somewhat 
larger  than  CC^' 

The  velocity  of  a  particle  on  a  typical  profile  is.  from  Equation  (52bj 


^  1  - _ £ -  (79) 

dt  t  e;;p(-z) 

so  that  the  velocities  of  two  particles  characterized  by  the  same  time, 
i,e.  belonging  to  the  same  profile,  are  equal  He  also  note  that,  since 
the  expression  for  the  velocity  can  be  written  in  terms  of  t  only, 
particles  characterized  by  different  masses  will  have  the  same  limiting 
velocities  although  the  corresponding  altitudes  will  be  different  This 
result  may  easily  be  seen  by  noting  that  all  points  on  a  profile  have  the 
same  velocity.  This  velocity  is,  therefore,  obtainable  from  the  velocity 


! 


at  maximum  concentration  which  is  eiitirely  a  drift  velocity  and,  there 
fore  (see  Equation  (10)  is  independent  of  mass. 

If  we  write  Equation  (52b)  in  the  form 


z  =  -  -enc  fat  (80) 

we  see  that  a  plot  of  z  against  fnt  for  a  fixed  value  of  c  should  result 

in  astraight  line.  Since  the  constant  c  enters  as  an  additive  term, 

lines  characterized  by  r'ffcrcnt  values  of  a  will  be  parallel  to  each 

other  The  constant  c  takes  on  a  particularly  simple  form  if  the  maximum 

2 


value  of  rij  is  used  in  Equation  (52a) 


For  this  case. 


a  so  that  here 


t 

c  =  a  and  the  lines  in  the  zvs--*nt  plot  are  now  character  ized  by  different 
values  of  a.  This  anticipated  result  is  verified  in  Figure  5  where  we  have 
plotted  on  semi- log  paper  the  time-altitude  history  for  several  values  of 
a  of  the  maximum  concentration  points  of  the  profiles  in  Figures  1  through 


4, 


The  initial  distribution  chosen  by  Banister  and  Davis  for  their 
Investigation  was  a  step  function; for  our  case,  a  delta  function  was 
used  It  is  easy  to  show  that  the  foregoing  results  are  independent  of 
the  initial  distribution-  For  a  delta  function,  initial  distribution 
6(z  ~  z^)  the  corresponding  solution  is  the  Green's  function  G(z,z^t) 
Using  the  properties  of  the  Green's  function,  we  may  write  the  solution 
n^(z,t)  corresponding  to  any  arbitrary  initial  distribution  f(z)  in  the 
form 
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EquaLiuii  (5za)  we  note  that  for  large  t,  G  splits  into  a  product 

of  a  function  F  of  z  and  t  and  a  function  -s:  of  z  , 

o 

G(z,z  ,t)  -  F(z,t)  rp(z  )  (821 

o  r  o 

We  then  have  for  large  t. 


n^(z,t)  - 

^  ~  1 
n 


F(z,t)  ^(z  )  f(z  )  dz 
o  o  o 


n^(z,t) 


=  A  F(z,t) 


where 


A  =  ,,(z  f(z  )  dz  (83b) 

o  o  o 

depends  only  on  the  parameters  H  and  a.  Thus,  the  solution  retains  the 
same  analytic  structure  in  the  region  of  interest,  independent  of  the 


initial  distribution. 
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